home *** CD-ROM | disk | FTP | other *** search
/ Fritz: All Fritz / All Fritz.zip / All Fritz / FILES / PROGMISC / PCSSP.LZH / PC-SSP.ZIP / MATINSL.ZIP / SIMQ.FOR < prev    next >
Text File  |  1985-11-29  |  3KB  |  123 lines
  1. C
  2. C     ..................................................................
  3. C
  4. C        SUBROUTINE SIMQ
  5. C
  6. C        PURPOSE
  7. C           OBTAIN SOLUTION OF A SET OF SIMULTANEOUS LINEAR EQUATIONS,
  8. C           AX=B
  9. C
  10. C        USAGE
  11. C           CALL SIMQ(A,B,N,KS)
  12. C
  13. C        DESCRIPTION OF PARAMETERS
  14. C           A - MATRIX OF COEFFICIENTS STORED COLUMNWISE.  THESE ARE
  15. C               DESTROYED IN THE COMPUTATION.  THE SIZE OF MATRIX A IS
  16. C               N BY N.
  17. C           B - VECTOR OF ORIGINAL CONSTANTS (LENGTH N). THESE ARE
  18. C               REPLACED BY FINAL SOLUTION VALUES, VECTOR X.
  19. C           N - NUMBER OF EQUATIONS AND VARIABLES. N MUST BE .GT. ONE.
  20. C           KS - OUTPUT DIGIT
  21. C                0 FOR A NORMAL SOLUTION
  22. C                1 FOR A SINGULAR SET OF EQUATIONS
  23. C
  24. C        REMARKS
  25. C           MATRIX A MUST BE GENERAL.
  26. C           IF MATRIX IS SINGULAR , SOLUTION VALUES ARE MEANINGLESS.
  27. C           AN ALTERNATIVE SOLUTION MAY BE OBTAINED BY USING MATRIX
  28. C           INVERSION (MINV) AND MATRIX PRODUCT (GMPRD).
  29. C
  30. C        SUBROUTINES AND FUNCTION SUBPROGRAMS REQUIRED
  31. C           NONE
  32. C
  33. C        METHOD
  34. C           METHOD OF SOLUTION IS BY ELIMINATION USING LARGEST PIVOTAL
  35. C           DIVISOR. EACH STAGE OF ELIMINATION CONSISTS OF INTERCHANGING
  36. C           ROWS WHEN NECESSARY TO AVOID DIVISION BY ZERO OR SMALL
  37. C           ELEMENTS.
  38. C           THE FORWARD SOLUTION TO OBTAIN VARIABLE N IS DONE IN
  39. C           N STAGES. THE BACK SOLUTION FOR THE OTHER VARIABLES IS
  40. C           CALCULATED BY SUCCESSIVE SUBSTITUTIONS. FINAL SOLUTION
  41. C           VALUES ARE DEVELOPED IN VECTOR B, WITH VARIABLE 1 IN B(1),
  42. C           VARIABLE 2 IN B(2),........, VARIABLE N IN B(N).
  43. C           IF NO PIVOT CAN BE FOUND EXCEEDING A TOLERANCE OF 0.0,
  44. C           THE MATRIX IS CONSIDERED SINGULAR AND KS IS SET TO 1. THIS
  45. C           TOLERANCE CAN BE MODIFIED BY REPLACING THE FIRST STATEMENT.
  46. C
  47. C     ..................................................................
  48. C
  49.       SUBROUTINE SIMQ(A,B,N,KS)
  50.       DIMENSION A(1),B(1)
  51. C
  52. C        FORWARD SOLUTION
  53. C
  54.       TOL=0.0
  55.       KS=0
  56.       JJ=-N
  57.       DO 65 J=1,N
  58.       JY=J+1
  59.       JJ=JJ+N+1
  60.       BIGA=0
  61.       IT=JJ-J
  62.       DO 30 I=J,N
  63. C
  64. C        SEARCH FOR MAXIMUM COEFFICIENT IN COLUMN
  65. C
  66.       IJ=IT+I
  67.       IF(ABS(BIGA)-ABS(A(IJ))) 20,30,30
  68.    20 BIGA=A(IJ)
  69.       IMAX=I
  70.    30 CONTINUE
  71. C
  72. C        TEST FOR PIVOT LESS THAN TOLERANCE (SINGULAR MATRIX)
  73. C
  74.       IF(ABS(BIGA)-TOL) 35,35,40
  75.    35 KS=1
  76.       RETURN
  77. C
  78. C        INTERCHANGE ROWS IF NECESSARY
  79. C
  80.    40 I1=J+N*(J-2)
  81.       IT=IMAX-J
  82.       DO 50 K=J,N
  83.       I1=I1+N
  84.       I2=I1+IT
  85.       SAVE=A(I1)
  86.       A(I1)=A(I2)
  87.       A(I2)=SAVE
  88. C
  89. C        DIVIDE EQUATION BY LEADING COEFFICIENT
  90. C
  91.    50 A(I1)=A(I1)/BIGA
  92.       SAVE=B(IMAX)
  93.       B(IMAX)=B(J)
  94.       B(J)=SAVE/BIGA
  95. C
  96. C        ELIMINATE NEXT VARIABLE
  97. C
  98.       IF(J-N) 55,70,55
  99.    55 IQS=N*(J-1)
  100.       DO 65 IX=JY,N
  101.       IXJ=IQS+IX
  102.       IT=J-IX
  103.       DO 60 JX=JY,N
  104.       IXJX=N*(JX-1)+IX
  105.       JJX=IXJX+IT
  106.    60 A(IXJX)=A(IXJX)-(A(IXJ)*A(JJX))
  107.    65 B(IX)=B(IX)-(B(J)*A(IXJ))
  108. C
  109. C        BACK SOLUTION
  110. C
  111.    70 NY=N-1
  112.       IT=N*N
  113.       DO 80 J=1,NY
  114.       IA=IT-J
  115.       IB=N-J
  116.       IC=N
  117.       DO 80 K=1,J
  118.       B(IB)=B(IB)-A(IA)*B(IC)
  119.       IA=IA-N
  120.    80 IC=IC-1
  121.       RETURN
  122.       END
  123.